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ABSTRACT 

The plasma density enhancements recently observed by the Large-Angle Spectrometric Corona- 
graph (LASCO) instrument onboard the Solar and Heliospheric Observatory (SOHO) spacecraft have 
sparked considerable interest. In our previous theoretical study of the formation and initial motion 
of these density enhancements it is found that beyond the helmet cusp of a coronal streamer the 
magnetized wake configuration is resistively unstable, that a traveling magnetic island develops at 
the center of the streamer, and that density enhancements occur within the magnetic islands. As 
the massive magnetic island travels outward, both its speed and width increase. The island passively 
traces the acceleration of the inner part of the wake. 

In the present paper a few spherical geometry effects are included, taking into account either the 
radial divergence of the magnetic field lines and the average expansion suffered by a parcel of plasma 
propagating outward, using the Expanding Box Model (EBM), and the diamagnetic force due to the 
overall magnetic field radial gradients, the so-called melon-seed force. It is found that the values of 
the acceleration and density contrasts can be in good agreement with LASCO observations, provided 
the spherical divergence of the magnetic lines starts beyond a critical distance from the Sun and the 
initial stage of the formation and acceleration of the plasmoid is due to the cartesian evolution of MHD 
instabilities. This result provides a constraint on the topology of the magnetic field in the coronal 
streamer. 

Subject headings: MHD — sun: corona — sun: magnetic fields — solar wind 



1. INTRODUCTION 

One of the most interesting findings of the LASCO 
instrument onboard the Solar and Heliospheric Obser- 
vatory (SOHO) spacecraft has been the observation of a 
continuous outflow of material in the solar streamer belt. 
An an alysis performed usin g a difference i mage tech- 
nique (jSheelev et al. I (jl997ft . iWang et all ((1998)) has 
revealed the presence of plasma density enhancements, 
called "blobs", accelerating away from the Sun. These 
plasmoids are seen to originate just beyond the cusps of 
helmet streamers as radially elongated structures a few 
percent denser than the surrounding plasma sheet, of ap- 
proximately IRq in length and O.li?© in width. They are 
observed to accelerate radially outward maintaining con- 
stant angular spans at a nearly constant acceleration up 
to the velocity of 200-450 kms^ 1 , in the spatial region 
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between about 5 and 30i?© . It has been inferred that the 
blobs are "tracers" of the sl ow wind, being carr ied out 
by the ambient plasma flow (Sh eelev et al. II1997I ). 

The solar streamer belt is a structure consisting of 
a magnetic configuration centered on the current sheet, 
which extends above the cusp of a coronal helmet. The 
region underlying the cusp is made up of closed mag- 
netic structures, with the cusp representing the point 
where separatrices between closed and open field lines 
intersect. Further from the Sun, at solar minimum, the 
streamer belt around the equator appears as a laminar 
configuration consisting of a thick plasma sheet with a 
density about 1 order of magnitude higher than the sur- 
rounding plasma, in which much narrower and complex 
structures are embedded. As a first approximation, mov- 
ing from the center of the streamer in polar directions at 
radii greater than the radius of the cusp, the radial com- 
ponent of the magnetic field increases from zero, having 
opposite values on the two sides of the current sheet. 
As far as the flow distribution is concerned, the fast so- 
lar wind originates from the unipolar regions outside the 
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Region of Interest 
(moving with the "blob") 




Expanding Computational Box 



Figure 1. Expanding Box Model: coordinate system and notations. The spatial coordinate aligned with the mean flow (streamwise 
direction) is denoted by y, x is the spatial coordinate along which the mean flow varies (cross-stream direction) and z is the coordinate 
along which physical fields are invariant (spanwise direction). We approximated the region of interest by a parallelepiped whose lengths 
are L along y, a and b along x and z respectively (not drawn in figure). The distance from the center of the Sun is R(t), which changes 
with time when the plasmoid moves outward. 



streamer belt, while the slowest flows are located at the 
center of the sheet. 

Explaining how both the slow and the fast component 
of the solar wind are accelerated is one of the outstand- 
ing problem in solar physics. Although the association 
between the slow sol ar wind and the streamer belt (e.g. 
iGosling et al. I ()1981[ )) and between the fast wind and the 
polar coronal holes is broadly recognized, the mechanism 
which leads to such acceler ations is still a matter of de- 
bate. lEinaudi et al. I (|1999t ) developed a magnetohydro- 
dy namic model (later ex tended to the compressible case 
in Einaud fet al. I ()2001[ )) that accounts for many of the 
typical features observed in the slow component of the so- 
lar wind: the region above the cusp of a helmet streamer 
is modelled as a current sheet embed ded in a broader 
wake fl ow. This and pr evious studies (jDahlburg et al. I 
(|1998() . lDahlburg| (|1998fl ) show that reconnection of the 
magnetic field occurs at the current sheet and that in the 
non-linear regime, when the equilibrium magnetic field is 
substantially modified, a Kelvin-Helmoltz instability de- 
velops, leading to the acceleration of density enhanced 
magnetic islands (Ein audi et al. Il200ll ). 



The problem we study has a roughly spherical sym- 
metry, but in our simulations we approximate the region 
beyond the cusp of a helmet streamer with a rectangu- 
lar box, and the geometry of the fields with a cartesian 
geometry, instead of a spherical one. To perform the 
numerical calculations we use a 2D numerical code with 
periodic boundary conditions in the streamwise direction 
(the spatial direction aligned with the mean flow) and 
non-reflecting boundary conditions in the cross-stream 
direction (the spatial direction along which the initial 
mean flow varies). A key feature, at the same time an 
advantage and a limitation of our analysis, is that thanks 
to the periodic boundary conditions that we use in the 
direction of the flow we are allowed to proceed in all our 
numerical simulations in a pseudo-Lagrangian fashion: 
i.e. the size of the computional box is fixed and then the 
time evolution of the spatial region moving with the box 
is followed in time, simply because what goes out from 
one side comes in from the other. In this case we follow 
the development of the magnetic island during its accel- 
eration outward along the radial direction. With respect 
to spatially developing studies of the solar wind our time 
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developing method is able to reach a higher resolution. 
In our case in fact the computational domain is smaller 
and we are allowed to increase the grid resolution along 
the current sheet in order to resolve the small scales de- 
veloping during the non-linear stage, and hence enabling 
us to study the development of a resistive-like instability 
that would be lost in a spatially developing study with 
a much lower resolution. On the other hand, so far no 
one has been able to study this phenomenon in spherical 
geometry, and we are planning to do it as the next step 
in the development of our model. 

In the present paper we introduce two major effects re- 
lated to the spherical symmetry of the problem we study. 
In what follows we suppose that when the magnetic is- 
land, the "blob", moves outward along the radial direc- 
tion it suffers a spherical expansion due to the approxi- 
mate spherical topology of the physical fields (Figured]). 
Although the subsequent algebra is quite complicated, as 
this approach enables us to attach the problem thanks 
to the availability of enough grid points, we still want 
to use the previous pseudo-Lagrangian time developing 
technique. This necessarily requires the use of a rect- 
angular computational box and fields with a Cartesian 
geometry, but we can include the expansion suffered by 
a parcel of plasma propagating outward by making the 
box expand at the correct rate. This is done by inserting 
a suitable force term in the momentum equation, a gen- 
eralization to the variable vel ocity regime of the Expand - 
ing Box Model develo ped bv iGrappin and Velli I $1996), 
IGrappin et al. I (|1993f ). as described in the following sec- 
tions. The other spherical effect included in the present 
paper is the so-called "melon-seed" force. This is a dia- 
magnetic force that a region of closed magnetic field lines 
suffers when it is embedd ed in a radially divergent back - 
gro und magnetic field (s ec Parker] (|1954f ). lParker I (|1957f ) 
and lSchmidt and Cargill I (120000 1. While in a spherical 
geometry this force would rise naturally, in our Cartesian 
approximation we take it into account by inserting a suit- 
able force term in the momentum equation, in analogy 
with the Expanding Box Model. 

Section [2] describes the governing equations as well as 
the initial and boundary conditions, the last ones de- 
scribed with more details in the Appendix. In section [3] 
we detail our numerical results. Section U contains a 
discussion of our results, as well as some concluding re- 
marks. 

2. PHYSICAL MODEL 

2.1 . Governing Equations 

In our analysis and simulations we use the compress- 
ible, isothermal, dissipative, three dimensional magneto- 
hydrodynamic (MHD) equations, that we write here in 
dimensionless form: 



dp „ 



(p«) = 0, 



(1) 



' du _ 

P I TTT + U ■ VU 



dt 



V P 



B 2 



B ■ VB (2) 

+ iv.£, (3) 

f) R 1 

_ = Vx(«xII) + — V"*, (4) 



supplemented by the equation of state P — pT. In the 
above equations p(x, t) is the mass density, u(x, t) is the 
flow velocity, P(x,t) is the thermal pressure, B(x,t) is 
the magnetic induction field, T(x,t) is the plasma tem- 
perature and £y = J^- + — §<%V • u is the viscous 

stress tensor. To render non-dimensional the equations 
we used the characteristic quantities L* , u* , p* , and the 
related quantities t* , T*, B* 



t* 



L* 



B* 



u*^f, (5) 



where m p is the proton mass. We use a simplified dif- 
fusion model, in which both the magnetic resistivity (rj) 
and the shear viscosity (p) are constant and uniform. A 
Stokes relationship is assumed, so that for the bulk vis- 
cosity A = (|) fjb. The kinetic and magnetic Reynolds 
numbers 1Z and lZ m are then given by: 



K = 



*u*L* 
P 



7? — 



Anu*L* 



rjc z 



(6) 



To ensure the solenoidality of the magnetic field, and 
supposing that there is no variation of the fields along z, 
we introduce the magnetic potential <f> defined by: 



B = Vx(^e 2 ) + B 2 
and replace equation {4]) with: 



f 



dt ' 
dB z 
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d 2 B z 
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du z 
V dy 

d 2 B z 
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(7) 

(8) 
(9) 
,(10) 



The equations solved in our numerical simulations de- 
pend on the problem that is studied: equations dU)- d3j) , 
((9^)- (ITU)) are solved when we do not take into account the 
expansion and the diamagnetic force, whereas when we 
consider these last effects some new force terms are added 
to the equations ([I])-©, (|§))- (ITU1) and a coordinate trans- 
formation is implemented for a matter of convenience, as 
described forward in this section. In any case the nu- 
merical problem is simplified by assuming that there is 
no variation in one of the spatial coordinates (z). In the 
following we refer to the spatial coordinate aligned with 
the mean flow as the streamwise direction (y) , the spatial 
coordinate along which the mean flow varies as the cross- 
stream direction (x), and to the remaining (spanwise) 
direction as z. The system has periodic boundary condi- 
tions in the streamwise direction, along which a Fourier 
pseudospectral method is used in the numerical compu- 
tations, and nonreflecting boundary condition s, achieved 
via the method of project e d characteristics (iThompson I 
(|1987l ). [Thompson 1 1)19900 . 1 Vanaiakshi et al. I (|1989l) h in 
the cross-stream direction, where a co mpact finite differ- 
ence scheme of the sixth order is used (jLele l fl992). This 
method, coupled with a hyperbolic tangent mesh stretch- 
ing around the current sheet, allows reasonable Reynolds 
numbers to be achieved without increasing the number 
of grid points dramatically. We solve the equations in a 
2D box whose dimensions are l x x l y = ±10.63 x 2ir/a, 
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where a is the streamwise wavenumber (the value of 
the box length along the x direction is the result of 
the mesh stretching). Time is discretized with a third- 
order Runge-Kutta method. In the simulations that we 
present in this paper we used a numerical grid with 
n x x n y — 401 x 256 points, and for all simulations the 
Reynolds numbers are 1Z = lZ m = 200. For a better 
des cription of the numer ical techniques used in our code 
sec lEinaudi et aT~l ((2001 . 

The region where the acceleration of the solar wind oc- 
curs would be more aptly described in spherical coordi- 
nates, but in order to perform our numerical simulations 
we do a Cartesian approximation, i.e. we approximate 
the fields topology with a Cartesian one. In this way 
we can use our pseudo-Lagrangian technique using a 2D 
code. As previously said we include the spherical expan- 
sion suffered by a parcel of plasma propagating outward 
by making expand our computational box at the correct 
rate. The region we study (see Figure [TJ) is a spherical 
sector of angular extent a in the cross-stream direction 
(x) and (3 in the span-wise direction, whose length is L 
along the stream-wise direction (y) and whose distance 
from the center of the Sun is R(t), where R changes 
in time while the plasma flows outward along the ra- 
dial direction. We approximate this region with a par- 
allelepiped whose cross lengths are respectively a and b 
along the directions x and z, and whose length in the 
radial direction is L. We suppose that when the plasma 
present in this region moves outward it suffers an aver- 
age expansion that we approximate with the expansion of 
the box (see Figure [I}. When following the evolution in 
the comoving frame of reference a non-inertial force term 
arises. In the Appendix it is shown that the force field for 
unity of mass that inserted in the momentum equation 
yields the correct expansion of the computational box is 
the following: 



f exp( X i t) 



R(t) 

R(t) 



(x,0,z). 



(11) 



In the previous equation the velocity of the box and its 
distance from the center of the sun are defined as: 



R(t) 



y Jo 



u y (0,y,0,t)dy, 



R(t) = R + R(t')dt' 



(12) 



(13) 



i.e. we define the velocity of the box as the average ve- 
locity in the y direction along the center of the computa- 
tional box. T he basic numerical algo r ithm was developed 
and te sted by IGrappin et al. I (|1993| ) , IGrappin and Veflil 
(1996) who showed that the EBM gives a correct descrip- 
tion of the evolution of waves (exact wave-action conser- 
vation) and turbulence in an expanding medium as well 
as capturing basic features of the expanding solar wind 
such as the formation of corotating interaction regions 
and forward and reverse shocks between high speed and 
low speed solar wind. The EBM does not conserve en- 
ergy because in an expanding flow work is done by the 
plasma. Angular momentu m and magnetic flux are c on- 
served exactly however (see lGrappin and Velli I ()1996l ) for 
a detailed description and tests of the physical model). 
The force term (|11|) introduces a variation along the x 



and z coordinates, but this explicit coordinate depen- 
dence cancels when we make a coordinate transforma- 
tion to a self-similarly expanding transverse coordinate 
x' (x, t) defined by 



Ro 

Wr 



y = y, 



Rq 

mt)' 



(14) 



It can be easily verified that in this way the expand- 
ing computational domain coordinate is transformed to 
a fixed computational grid. The equations solved by our 
2D numerical code and all the mathematical details are 
extensively explained in the Appendix. 

When we make our Cartesian approximation we ap- 
proximate the radially diverging magnetic field lines with 
straight ones, losing the radial nonhomogeneity of the 
solar magnetic field. An effect that is lost with this 
approximation is the so-called "melon seed" force, a 
diamagnetic force acting on a r egion of closed mag- 
netic fie ld lines. It is known ( s ee ISchmidt and Cargill I 
(20001 . IPneuman and Careilil (|1985h . IParker I (|1954T > 
and lParker I (]1957l )) that a region of closed magnetic field 
lines embedded in a background radially diverging mag- 
netic field is subject to such force. Modeling this region 
as a prolate spheroid of volume V elongated in the radial 
direction it can be demonstrated that the value of the 
total force acting on the spheroid is approximately: 



Fm =~V 



\- 


(BIX 


dR 





(15) 



where B e = B c q (^) is the background magnetic field, 
-B e o being constant. In our model axial invariance is as- 
sumed, i.e. there is no variation in the z direction in the 
Cartesian approximation. Hence the prolate spheroid is 
modeled as a cylinder with an elliptical cross section, and 
the same force for unit volume is used. The force given in 
equation (|15[) is the total force acting on this plasmoid. 
The force per unit volume derived from this equation and 
used in our computations is given by: 



f m {x,y,t) = Q(x,y,t) 



B e0 

2vri? n 



Ro 
R 



that written in non dimensional form is given by: 



2 #e0 f Ro^ 5 



(16) 



f m ( x , y ,t) = e(x,y,t)^\^j e y , (17) 

in which is a field whose value is — 1 inside magnetic 
island and = outside. In this way the field (| 
integrated over the volume of the plasmoid gives as a 
result the total force (fT5j). 

The region above the cusp of a helmet streamer is as- 
sumed to be in local hydrodynamic equilibrium, there- 
fore gravity, magnetic curvature forces and the flow all 
balance out. Once density fluctuations arise, the gravita- 
tional force reappears as an Archimedes buoyancy force 
acting on regions with a different density compared to 
the average. This force field can be written in non- 
dimensional form as 



fa -Lg^ 



where ka = 



GM Q 

l* (u*y 



(18) 



G being the Gravitational constant and Mq the solar 
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mass. Taking as characteristic length and velocity, the 
width of the wake L* = ay ~ Q.IRq and the velocity of 
the fast solar wind u* — ~ 6 • 10 5 ms" 1 , we obtain 
the value kc — 5.3. We can now compare the relative 
strength of this buoyancy force with that of the melon- 
seed force (|17[) . As shown in the next paragraph A ~ 2.5, 
Ro = 60 (i.e. 6 Rq in conventional units, and supposing 
A/? = 1, we obtain near the cusp of the helmet streamer 
(R = Ro): 



frn 



Ml 

Ro 

G R-2 



100. 



(19) 



i.e. the buoyancy force is negligible with respect to the 
melon seed force near the cusp of the helmet streamer 
and in the region of interest of the simulations presented 
in the following paragraphs. 

In the next section we introduce as initial conditions 
the fields (|22 p -(|24 |) . These fields are a one-dimensional 
equilibrium for the region beyond the cusp of a helmet 
streamer, i.e. the variation of the fields is only along 
the cross-stream direction (the x coordinate) and So- 
lar Corona Reynolds numbers are so high that they are 
substantially an equilibrium solution of equations (fT)- 
(|9j) over the time scales of our simulations. But in our 
computational domain Reynolds numbers are limited by 
resolution to values around 1Z = lZ m = 200, values which 
are so small that unphysical diffusion of equilibrium fields 
would take place. In the foregoing equations the pattern 
of the diffusion terms is 



df 

at 



n dx 2 ' 



(20) 



where / is a generic field. As the variation of the equi- 
librium fields is only along the x coordinate, i.e. as 
a Fourier series in the y direction they have only a 
k = mode, we avoid such an unphysical diffusion using 
instead of eq. (|20p the following one: 



df 

at 



i a 2 



f-fo 



(21) 



In this way the modes k — of the Fourier series in 
the y direction, that we indicate with /o, do not diffuse. 
This could have potentially given rise to development 
of excessively small scales in x on the k = harmonic 
but empirically the numerical filtering used along this 
direction has been enough. 

2.2. Initial Conditions 

We model a planar section of the solar streamer belt, 
in the region beyond the cusp of a helmet streamer, as 
a magnetic current sheet of thickness embedded in 
the center of a broader wake flow of thickness ay. The 
plasma flows at the speed of the fast solar wind at the 
edges and at a much lower velocity at the sheet. Taking 
as characteristic length the thickness of the wake L* = 
ay, and as characteristic density and velocity those of the 



fast wind p* 



the basic fields, written in 



non-dimensional form, are then given by: 
uo y (x) = 1 — scch (x) , 



(22) 



R=200 



0): 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 

a 



Figure 2. Dispersion relation for the resistive mode with A = 1.5, 
R = 200. The growth rate (a;;) as a function of wavenumber (a) 
is displayed (in non-dimensional units). 



P =1 = T = T1F>. 



1 

M 2 ' 



(24) 



where 8 = — is the ratio of the two widths, and A and 

a B ' 

M are respectively the Alfven number, the ratio between 
the Alfven speed (caoo) an d the flow speed (moo), and the 
sonic Mach number, the ratio between the flow speed 
(ttoo) an d the sound speed (c soo ), of the system at the 
edges, i.e. of the fast wind: 



.4 



CAc 



OO^C 



M = 



— 2£_. (25) 



Bo y (x) = Atanh (Sx) 



B 0z (x) = A sech(fa) , (23) 



Our wake-current sheet model is then characterized by 
three parameters, A, M, and S, which vary consider- 
ably in the solar corona as a function of the distance 
from the Sun. In this paper we are interested in study- 
ing the time and spatial evolution of a parcel of plasma 
starting from the region immediately beyond the cusp 
of a helmet streamer. In this region the typical Alfven 
velocity c^oo is approximately 1000-1500 km s" 1 , the 
sound speed c soo is about 100 km s^ 1 , and the velocity of 
the fast solar wind Uoo is ~ 600 kms^ 1 . It follows that 
A = c Aoo / Uoo - 1500/600 = 2.5, and M = Moo/c soo ~ 
600/100 = 6. The plasma j3 is also a useful quantity, 
defined as /3 = 2p/\B\ 2 = 2c sao 2 /c Aoo 2 = 2/ (AM) 2 . 
In this case we have /3 ~ 0.0089. In our numerical 
calculations we use values of these parameters that are 
close to these estimates but computationally more acces- 
sible. Thus, in our runs A — 1.5 and M — 2 (hence 
f3 = 2/9 = 0.22). Unfortunately there is no direct ob- 
servation of the heliospheric current sheet close to the 
Sun (r < lOi?©), hence there is no measurement of the 
parameter S. But from the topology of both the mag- 
netic and the velocity fields it is reasonable to suppose 
S » 1, i.e. the width of the wake flow much big- 
;er th an that of the current sheet (see iStribling et al. I 
1996fh lEinaudi et all (|1999| )). Such a high value of S 
is not computa tionally accessible and on the other hand 
we have seen in lEinaudi et al. I (|1999f) that the impact of 
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Figure 3. Run A — No expansion or melon-seed force is included in this simulation. Time evolution of magnetic islands is showed. 
Magnetic field lines (white lines) and electron density filled contours (colours show the value of the density field in accordance with the 
colour scale) at selected times show density enhanced magnetic islands forming and accelerating outward. The x and z axis are respectively 
along the vertical and horizontal direction and are measured in solar radii units (Rq). ni, is the average electron density at the top and 
bottom z boundaries (x = ±1.06 Rq), it doesn't change significantly because the box is not expanding. A ~ 5 hours time should be added 
to the times showed to keep count of the rescaled linear time-scale. 




. 

-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 



Figure 4. Run A — Streamwise component of the velocity (u y ) 
averaged in the y direction as a function of x (the cross-stream 
direction) at selected times. When the nonlinear regime develops 
(after t ~ 58 m), while magnetic reconnection continues to occur, a 
Kelvin-Helmholtz-like instability is triggered, leading to the accel- 
eration of the central part of the wake up to the value ~ 500fcm s _1 . 



the value of 8 is mainly to change the critical value of 
A below which the dynamics is dominated by the flow. 
Therefore in this paper we have chosen 6 — 5, which 
corresponds to a magnetically dominated regime. 

3. RESULTS 

This section presents the results of our numerical sim- 
ulations. In £13.11 we show for reference the result of a 
numerical simulation in which both the melon seed force 
and the spherical expansion are neglected. In H3.2[ in 
order to show the influence of the expansion on the dy- 
namics of the wake-current sheet system, we present the 
results of a numerical simulation in which the geometri- 
cal expansion is implemented while the melon seed force 
is neglected. At last, in H3.3I we show the results of a 
simulation in which both the melon seed force and the 
spherical expansion are implemented. 

If there was no magnetic field (A — 0) the wake-current 
sheet system (equations ([!?!?)) - ([iM]) ) would be unstable for 
a Kelvin- Hclmholtz-like instability. Two unstable modes 
exist in this limit: a varicose (sausage-like) mode and a 
sinuous (kink-like) mode. The stabilizing effect of a mag- 
netic field parallel to the flow o ver the Kelvin-H e lmhol tz 
instability is well known (e.g. [Chandrasckh ar] (119611)') . 
and in fact ou r linea r studies (jDahlburg et al. I (|1998|) . 
lEinaudi et all (fl999)) have found that increasing the 
value of the Alfven number A beyond a critical value 
the two ideal modes are stabilized. On the other hand 
increasing the strength of the magnetic field leads the 
system to be unstable to a tearing-like instability. This 
resistive varicose mode has the same spatial symmetry 
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Figure 5. Run B — Expansion but no melon-seed force included in this simulation. Time evolution of magnetic islands is showed. 
Magnetic field lines (white lines) and electron density filled contours (colours show the value of the density in accordance with the colour 
scale) at selected times show density enhanced magnetic islands forming and accelerating outward. The x and z axis are respectively along 
the vertical and horizontal direction and are measured in solar radii units (Rq). The initial distance between the box and the center of 
the Sun (Ro) is 6Rq, during time evolution this distance increases, and at the selected times showed in the Figure the values of R(t) are 
given by: 6.02, 6.2, 7.01, 8.11, 9.32, 11.86 Rq. rif, is the average electron density at the top and bottom z boundaries, as the box expands its 
value decreases as (Ro/R) 2 . A ~ 5 hours time should be added to the times showed to keep count of the rescaled linear time-scale. 



of the varicose fluid mode and, while in the linear regime 
the dynamics is not substantially modified by the pres- 
ence of the flow, developing like a classic tearing insta- 
bility, in the non-linear regime the presence of the flow 
has a dramatic influence. In fact, as the perturbed mag- 
netic and velocity fields attain finite amplitude, the re- 
sultin g nonlinear stres ses along the magnetic island bor- 
ders |Dahlburg1 ll998f ) lead to a transfer of momentum 
between the flow and the magnetic island. As a result the 
magnetic island and the central part of the wake are ac- 
celerated. The coupling between magnetic reconnection 
and Kelvin-Helmholtz instability that dev elops in the dy- 
nami cs of this mode has been studied in (Ei naudi et al. I 
[19991) . 

In all calculations a small a mplitude perturb a tion o f 
the resistive varicose type (see fPahlbur g et ai~l fl998)) 
at the wavenumber a — 0.42 was added to the initial 
state. Although our numerical code was developed to 
study the non-linear dynamics, we can estimate (via a 
fitting during the linear regime) the growth rates of the 
perturbations at various wavenumbers. Figure [5] shows 
the dispersion relation for resistive instabilities obtained 
with this method. It is evident from this Figure that 
we have not chosen the fastest growing wavenumber for 
our nonlinear calculations. The reason is related to the 
coupling of the magnetic reconnection instability and of 
the Kelvin-Helmholtz instability. The initiation of the 
Kelvin-Helmholtz instability occurs when the resistive 
instability is well into its nonlinear regime. The maxi- 



mum dej^edvajue imcy^ well k nown fact 

(e.g. iRavl (jl982fl . lFerrari and Trussoni I (|1983| )) that the 
Kelvin-Helmholtz instability is stabilized for wavelengths 
approaching the thickness of the shear layer (in our case 
the amplitude of the wake flow ay), in this way the wave- 
lengths that give rise to a more efficient acceleration of 
the wake are the longer ones (compared to the ampli- 
tude of the wake). The minimum desired value for a 
is determined by the number of subharmonics we want 
to include, to allow for vortex merging. For this rea- 
son, based on the Kelvin-Helmholtz dynamics, we have 
picked up the value a = 0.42 and not the wavenumber of 
the fastest growing resistive mode. A most comprehen- 
sive study of this mechanism is underway, in particular 
a multiwavelength one. 

While the numerical computations have been carried 
out using non-dimensional equations, for a matter of clar- 
ity the results will be presented in the ordinary dimen- 
sional units. As in this paper we are interested in study- 
ing the time and[ spatial evolution of a parcel of plasma 
starting from the region immediately beyond the cusp 
of a helmet streamer, i.e. 4 — 6 Rq, we have chosen 
as characteristic quantities those characterizing this re- 
gion. We have chosen as characteristic length the thick- 
ness of the wake L* = O.li?©, as characteristic velocity 
that of the fast wind u* = 600 kms" 1 , and as char- 
acteristic density the current shee t electron density at 
~_5i?0, i-e. n = 5 • 10 5 cm~ 3 (|Guhathakurta et al. I 
1996). From these quantities we derive the characteristic 
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Figure 6. Run B — Strcamwise component of the velocity (u y ) 
averaged in the y direction as a function of x (the cross-stream 
direction) at selected times. Expansion effects imply a weakening 
of cross-stream gradients of magnetic and velocity fields, leading 
to a lower acceleration profile. 




Figure 7. Run B — Average streamwise component of the ve- 
locity (< u y >) at selected x n values (x n representing the non- 
expanding numerical grid related to the expanding physical grid x 
by the relation x = x n ■ R(t)/Ro) as a function of R, the distance 
of the expanding box from the center of the Sun. 



time t* = L*/u* ~ 116 s. A particular caution should 
be taken when considering the times presented in the 
following simulations. Our numerical simulations have 
been in fact performed using Reynolds numbers values 
(R = 200) much lower than coronal values. While most 
of the nonlinear dynamics is due to ideal processes (i.e. 
Kelvin-Hclmholtz instability) which are little affected by 
dissipation, the growth rate of the tearing-like instabil- 
ity developing in the linear regime is strongly affected by 
the Lundquist number, hence the time-scale of the lin- 
ear stage must be rescaled. With a Lundquist number of 
the order S = 1 9 the linear phase s hould last approxi- 
mately 5 hours (jEinaudi et al. I I1999). This time should 
be added to the times of the simulations presented in the 
following paragraphs. 

3.1. Run A — No expansion or melon- seed force 

In this simulation we consider neither the expansion 
nor the diamagnetic force. While the parameters are 
somewhat different, we repr oduce the esse n tial e volu- 
tionary features repo rted by lEinaudi et al. I ()1999D and 
lEinaudi et al. I (|2001l ). Figure |3] shows the formation and 
acceleration of density enhanced magnetic islands. These 
magnetic islands form with streamwise length equal to 
the perturbation wavelength. The magnetic island cross- 
stream width grows from a small amplitude in the linear 
regime to the order of as = 0.02 i? Q , the width of the 
current sheet, at the beginning of the non- linear stage 
(t ~ 77m). While magnetic reconnection continues to 
occur in the nonlinear regime, the Kelvin-Helmholtz in- 
stability also is triggered, leading to the acceleration of 
the central part of the wake (Figure |4j up to the value 
~ 500 kms^ 1 at t = 309m 

The value of the density enhancement {n max /nb de- 
fined as the ratio of the maximum density inside the mag- 
netic island over the background density value, steadily 
increases from an initial uniform value (n max = n = 
5 • 10 5 cm~ 3 ) up to a value which is roughly 20 times 



greater (n max ~ 9.7- 10 6 ) at t = 309m. Whereas the for- 
mation and evolution of this density enhanced magnetic 
island is very similar to that of the observed "blobs" 
(i.e. we have a dense, accelerated plasmoid), both the 
density enhancement and the acceleratio n profile result 
in lar g er values than the o nes observed (|Sheelev et al. I 
(|1997D . IWang et al.1 (|1998f )). This has been one of the 
main reasons to study the more realistic configurations 
described in the next two sections. We will show in the 
next section that when expansion effects are included 
such a peaked density does not occur because, as ex- 
pected, while the magnetic island travels outward there 
is a rarefaction of the island due to the expansion. 

3.2. Run B — Expansion but no melon-seed force 

This run provides us with information on the effects of 
the geometrical expansion on the dynamics of the wake- 
current sheet model for slow solar wind formation and 
acceleration. In this section we neglect the diamagnetic 
force. We will account for this force in the next section. 

We used as starting radius i?o = 6i?©, approximately 
the radius of the region located near the cusp of a hel- 
met streamer. During the time evolution the radius R 
increases its value up to R ~ 12 R@ at t = 309m. From 
Figure [1] we can clearly notice that while the plasmoid 
moves outward the box suffers an expansion. Its volume 
increases as (R/Rq) 2 leading to an average decreasing of 
the mass density as (Rq/R) 2 , as the values of the back- 
ground density rib indicated in Figure [T] show. 

The dynamical evolution (Figure [SJ for the expansion 
case exhibits some crucial differences from that described 
in the previous section for the non-expanding case. Again 
we have the formation of a density enhanced magnetic 
island, but as now expansion is occurring the average 
density of matter decreases as (Rq/R) 2 , and the den- 
sity enhancement is lower. The density enhancement 
nmax/nb for the selected times shown in Figure [5] range 
between 2.5 and 5. This value is closer, with respect to 
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Figure 8. Run C — Expansion and melon-seed force implemented. Time evolution of magnetic islands is showed. Magnetic field lines 
(white lines) and electron density filled contours (colours show the value of the density in accordance with the colour scale) at selected times 
are showed. The value of the melon-seed force is very small during the linear stage, but as non-linear stage occurs, both Kelvin-Helmholtz 
and the diamagnetic force accelerate the island. The x and z axis are respectively along the vertical and horizontal direction and are 
measured in solar radii units (Rq). The initial distance between the box and the center of the Sun (Ro) is 6R@, during time evolution 
this distance increases, and at the selected times showed in the Figure the values of R(t) are given by: 6.58, 7.82, 9.46, 11.4, 13.5, 18.3 Rq. 
n b is the average electron density at the top and bottom z boundaries, as the box expands its value decreases as (Rq/R) 2 . A ~ 5 hours 
time should be added to the times showed to keep count of the rescaled linear time-scale. 



the non-expanding simulation presented in the previous 
paragraph, to the observed one, w hich is of the order o f 
7±2% of the back ground K corona (|Sheelev et al. lH99l . 
Dahlburg ! (|1998| ) showed that the transfer of energy to 
the perturbed fields from the background magnetic and 
velocity fields depends strongly on the cross-stream gra- 
dients of these fields. On the other hand the expansion 
of the box implies flow velocities orthogonal to the radial 
direc tion which lead to a weake ning of these field's gradi- 
ents (jGrappin and Velii"lll996[ ) and as a consequence to 
a lower acceleration profile, shown in Figure [6] 

In figure we plot the average streamwise velocity at 
selected x n values (here x n is the non-expanding numer- 
ical grid, related to the expanding physical grid x by the 
relation x = x n ■ R(t)/Ro) as a function of the distance 
between the box and the center of the Sun, which reaches 
the value R ~ 12i? Q in conventional units, at t ~ 309m. 
This velocity profile is in good quantitative agreement 
with the one obtaine d by observat i ons ( compare with 
Fig. 4 on page L166 of IWang et al~l (|1998H ). 



3.3. Run C — Expansion and melon- seed force 



In this section we present the results of a numerical 
simulation in which both the geometric expansion and 
the diamagnetic, or melon seed, force are implemented. 

As detailed at the end of $2. li the region above the cusp 
of a helmet streamer is assumed to be in local hydrody- 
namic equilibrium, therefore gravity, magnetic curvature 
forces and the flow all balance out. Once density fluc- 
tuations arise, the gravitational force reappears as an 
Archimedes buoyancy force acting on regions with a dif- 
ferent density compared to the average. We have shown 
that this buoyancy force is negligible with respect to the 
melon-seed force in the region of interest of our simula- 
tion. 

The expression for the melon seed force (jTTJ) holds for a 
magnetic island embedded in an external magnetic field 

B e = B e0 (-^-) . We compute the value of B e as the 
value of the magnetic field at the outer border of the 
magnetic island. Hence the value of the melon-seed force 
is very small during the linear stage, when the magnetic 
island is smaller than the current sheet width and the 
ambient magnetic field is small too. Figure [8] shows the 
time evolution of the system. During the linear regime 
the diamagnetic force is negligible, so that the dynam- 
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Figure 9. Run C — Streamwise component of the velocity (u y ) 
averaged in the y direction as a function of x (the cross-stream 
function) at selected times. The value of the melon-seed force is 
very small during the linear stage, but as non-linear stage occurs, 
both Kelvin-Helmholtz and the diamagnetic force syncrgistically 
act to accelerate the magnetic island leading to an over-acceleration 
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Figure 10. Run C — Average streamwise component of the ve- 
locity (< u v >) at selected x n values (x n representing the non- 
expanding numerical grid related to the expanding physical grid x 
by the relation X — Xn ■ R(t)/Ro) as a function of R, the distance 
of the expanding box from the center of the Sun. 



ical evolution is very similar to the ones presented in 
the previous simulations. But as the non-linear stage 
occurs, both the Kelvin-Hclmholtz and the melon-seed 
force, synergistically accelerate the island, leading to an 
over-acceleration of the plasmoid, as shown in Figures 

On the other hand, the expression of the diamagnetic 
force flTTl) holds for a plasmoid in a radial divergent back- 
ground magnetic field. Indeed we have implemented this 
diamagnetic force supposing that the magnetic field in 
the region just beyond the cusp of an helmet streamer 
diverges radially. Unfortunately there is no direct ob- 
servation of the topology of the magnetic field in this 
region. But as the resu l ting a c celeration profile is not 
observed (jSheelev et al. I (|1997t ). lWang et al. I (|1998D ) we 
conclude from the results of this simulation that in such 
a region the magnetic field line do not diverge (in other 
words, the melon seed must be much smaller than the 
semi-empirical estimate used in our simulations). Both 
more refined observations of this region and theoretical 
studies, including high resolution simulations in a fully 
spherical geometry, are desirable to establish what kind 
of magnetic topology we have in this region. 

4. DISCUSSION 

In this paper we have refined the model presented in 
lEinaudi et all (|1999f ). lEinaudi et al. I (|2001[ ). for the for- 
mation and acceleration of the slow component of the 
solar wind within a coronal streamer. In particular, we 
have been concerned with the study of the initial stage 
of the process. The main results of the present paper can 
be summarized as follows: 

1. The expanding box generalization of the wake- 
current sheet model exhibits the formation and out- 



ward acceleration of density enhanced magnetic is- 
lands. The plasmoid width, length and accelera- 
tion profiles are in good agreement with LASCO 
observations. 

2. The geometrical expansion by itself has a stabiliz- 
ing effect on the Kelvin-Helmholtz instability de- 
veloped in the non-linear regime. It also is respon- 
sible for a rarefaction of the accelerated plasmoid. 

3. From our simulations there is a strong indication 
that the topology of the magnetic field lines in 
the region close but beyond the cusp of an hel- 
met streamer is not radially divergent. This diver- 
gence would lead to as yet unobserved acceleration 
and density profiles due to the diamagnetic force. 
Note however that other effects (i.e. drag from the 
denser current sheets) might be important as well. 

In conclusion, these computations reinforce the sce- 
nario that the structure and acceleration of the slow solar 
wind can be attributed primarily to the intrinsic insta- 
bility and subsequent evolution of the underlying wake- 
current sheet system. To reproduce the observational 
details properly, further theoretical work and additional 
solar data are needed. From the theoretical point of view, 
we need a study of more complex and realistic basic con- 
figurations, where an initial density profile peaked on the 
current sheet is present, and a non force- free magnetic 
field is chosen. 
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APPENDIX 

As previously said we suppose that when the magnetic island moves outward the plasma present in the region in 
and around it suffers a spherical expansion due to the approximate spherical topology of the physical fields (Figured]). 
When a particle or a fluid element moves outward along the radial coordinate r at fixed azimuthal and polar angles, 
the ratios x/y and z/y are constant in time. For the sake of simplicity consider a system of spherical coordinates where 
9 is the polar angle from the y-axis (with < 9 < 7r), <j> is the azimuthal angle in the zx-plane from the z-axis (with 
< 4> < 27r), and r is the distance from the origin. As usual the relations between the spherical and the cartesian 
coordinates are given by: 

y = r cos 9, z = r sin 9 cos (f>, x — r sin 9 sin <f>. (1) 

It can be easily verified that the ratios x/y and z/y depend only by the angles 9 and <f>; hence for a particle or a fluid 
element moving along the radial coordinate these ratios are constant in time. In our case y = R(t) and for such a fluid 
element we can write: 

x(t) _ xo_ z(t) _ Zq_ 

R{t) ~ Ro ' R(t) ~ Rq' U 
where Xq, z and i?o are the initial values of these quantities. Differentiating with respect to time we have 

m - ^R® = m = = (3) 

This means that the cross-component of the velocity is given by: 

nx = M(x(t),Q,z(t))' (4) 

This holds for a particle or a fluid element moving outward along the radial coordinate. To apply this to our subject 
we replace this velocity with the fluid velocity and look for a force able to cause this fluid motion. In other words, in 
our picture a fluid that is expanding spherically is characterized by a velocity field whose cross-component is given by: 

u ± (x,t) = ^(x,0,z\ (5) 

Hence we introduce a new force in the momentum equation corresponding to the spherical expansion effects, i.e. the 
force term whose solution is given by the previous term. To find this force term we simply insert in the momentum 
equation the velocity field u±(x,t) which gives us the force field for unity of mass 

t) = ^ + (ul • V) uj. = ||| (x, 0, z) . (6) 

We have chosen to define the velocity of the box R as the average velocity in the y direction along the center of the 
computational box (see eq. fT2|) . and its distance from the center of the Sun R accordingly (see eg. fT3|) . 

In order to perform the numerical simulations the computational box must be fixed. For this reason we consider the 
coordinate transformation x' {x, t) defined by: 

X= W) X ' v=m Z =W) Z - () 

It can be easily verified that this coordinate transformations brings the physical dominium expanding with the law ^ 
into a fixed one. This coordinate transformations implies a transformations of the field. To the field /(a;, t) corresponds 
a new field g(x',t) defined by 

f(x,t)=g(x'(x,t)). (8) 

The following relations can be easily verified: 

^(x,t) = ^( x '( Xl t),t) - [U(x'(x,t),t) ■ V']g(x'(x,t),t) (9) 
V/(aB > f) = V / s(a:'(a:,i),t) (10) 
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U(x',t) = 



V 
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Rq 



-x',0 



Rq 



Rq 



(x',0,z>) 



( R d d R Q d 



\R{t) dx" dy" R{t) dz' y 
Defining the transformed fields p, T, u, B, with a primed index, like for example 

u(x,t) = u'(x'(x,t)), 
and noticing that the transformed expansion force field © has the following form 



f exp( X J 



R(t) (_j 

R(t) 



(x',Q,z')=^{x',t), 



(11) 
(12) 

(13) 

(14) 



the MHD equations (for clarity we explicitly include the isothermal equation) can be written in the new coordinate 
system and with the transformed fields as: 



| = ([/.v') P '-v'-(A') 



du' 

~dt 



([/• V'K +/'-(«'• V>' 



exp 

OB' 
~~dt 



VP' 



(B' ■ V) B' 



p'TZ 



(U • V) B' + V x («' x B') + V' 2 B' 



or 



(U • V) T' — (u' • V') T' 



where 



4- = - 3 v ' • <j = W + pl = p' r - 



In order to simplify the previous equation notice that 

R 



V ■ U = 2 



R 



(A ■ V) U = ^ & ± A , 



(15) 
(16) 

(17) 
(18) 

(19) 
(20) 



where A is a generic vector field and <^j_ is the projector in the cross-stream direction, i.e. £^j_A — (Aj.,0, A3). 
To render more simple our equations we introduce the new velocity field u(x',t) = u'(x',t) — U(x',t) the previous 
equations can be rewritten as 



du 
~dt 



r 



exp m J 



dU\ „, N . R _ 1 

■ — (it ■ V ) u — —0>±u 

R p' 



VP' 



B 



-(B'-V')B' 



— V-£' 

p'TZ q 



no/ D D 1 

i- = Vx(ixfl')- 2^B' + > X B' + ^V' 2 !?' 
Ci K it rt m 



dV 
~~dt 



(u-V')T' 



(21) 
(22) 

(23) 
(24) 



In particular the first term between round brackets in the momentum equation (|22[) cancels thanks to the equation (I14p . 
In this way we have a problem that numerically is invariant along the z direction and hence can be solved with a 2D 
numerical code. 

To simplify further the equations we introduce the average scalings of the physical fields obtained from the matter 
conservation law and the Alfven theorem. The new fields are the ones with the hat: 



Rc, 



R 



Rc 



R 



P = Kr A B x =[-1 B x , B y =\-^\ B y , B z = B z , 



Ri 



R 



P(; 



R 



(25) 



and for a matter of convenience we introduce the scalings 
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and 



v=|-,-,- 

v dx' ' dy' ' dz' 



Notice that for two functions / and g for which holds 



for the time derivative it follows 



Furthermore it can be easily shown that 



dg 
dt 



Ro 
~R 



df R 



(27) 



(28) 



(29) 



u ■ V = u ■ V , V • u = V • u (30) 

At last, the equations that we solve with our 2D numerical code are given by (in what follows all the fields components 
should have a hat, that we omit for the sake of simplicity): 



dp ( dp 



du x 
dt 



( du^ 
\ dt 



du x 



dt ' 

Rq 
R 



dt 



dp du y 



Uy dy P 'dy 



B y dB x 2R 1 
p dy R x plZ 



+ 



i d 2 



+ 



3 9a; 2 3 dxdy dy 2 



du y 
~dT 



dUy 

~dt 



fdT + Tdp\ B x dB x B z dB z I 



dUy 

v dy \dy ^ p dy J p dy p dy ^ pK 



d 2 u v 1 d 2 u x 



4d 2 u v 



du z 
dt 



( ^1 
V dt 



du z 

ly ~dJ 



Ro\ By dB z 



p dy 



2R 1 

-r Uz + 7r 



dx 2 3 dxdy 3 dy 2 
+ 



dB z 
dt 



dB z 
dt 

dBy 

dt 



dB x 
dt 



R du v 
dy 



dB,. 



B 



dt 



du z 
V dy 

dB 



dB z 



1 



V dy ^ n n 



dx 2 dy 2 
d 2 B z d 2 B 



dx 2 



dy 2 



v 



dT 
~dt 



V dy 

/dT 
\~dt 



d 2 B„ d 2 B„ 



d 1 

- { U X By-UyB X ) + — 



dx 2 
dT 

x d 2 B x 



dx 2 



+ 



dy 2 



d 2 B x 



dy 2 



where with the symbol ( ) x we indicate the non-dimensional terms containing only the partial x derivatives: 



du x 
~dt 



= -u x 



du x 
dx 



Ro 
~R 



( ® u y 
\~dt 

du z 
dt 

dB z 
dt 

fdBy^ 

\ dt 



dp) 




dp 


du x 


dt J 


= -Ux 

X 


dx 


^ dx 


dT 


+ Tdp\ 


( 




dx 


p dx J 





dB z 
dx 



+ 



Ro 
~R 



' dB y 
By -dx~ 



= -u 



dUy 

dx 

du z 
x dx 



-Ro \ B x dBy 



p dx 



. Ro\ B x dB 



R 

du 



p dx 
dB z 



dx dx dx 

du x du y dBy 

= —By— h B x — u x — — , 

dx ox dx 



(31) 
(32) 
(33) 
(34) 
(35) 
(36) 
(37) 
(38) 

(39) 
(40) 
(41) 
(42) 
(43) 
(44) 
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dT 

~at 



dT 

dx 



(45) 



In equations (|5T|) - (I551) also the diffusive terms would be affected by expansion factors, but we neglect these factors 
because in numerical codes diffusive terms play a numerically stabilizing role, which would be affected by such expansion 
factors. 

From the terms (|55|) - (f4"5"j) one constructs the charact eristic polynom i al for the x derivative terms, that we use 
to implement the nonreflecting boundary conditions (see lEinaudi et al~l (|2001)). This yields the projected entropy, 
Alfven, fast, slow characteristic speeds for the Expanding Box Model A in non dimensional form as 



X a = u x , \ a ±=u x ±£a' x , X f i 
where a' and c' are the modified Alfven and sound speeds 

/ -^o , Ro B x 

C= T VT > a * = T 4f> ' °* 

and £ = a' x /\a' x \. The fast and slow speeds are defined by: 



= u x ± f, X s ± =u x ±s' 



Ro 
R 



B v 



Ro -B* 

R VP ' 



(46) 



(47) 



,/2 _ 1 

1 _ 2 



(c' 2 + a' 2 ) + V / (c' 2 +a' 2 ) 2 -4<V 2 



(48) 



( C ' 2 + «' 2 )-\/( 



; ' 2 + a' 2 ) 2 -4a' 2 c' 2 



(49) 



Indicating with a[ , a' 2 and a' ± the quantities 



< f' 2 -< 2 
f' 2 - s' 2 



i 

a 2 



I f' 2 - c' 2 
f' 2 ~ s' 2 



a y 2 + a' z 2 



(50) 



the projected characteristics are then given by the following expressions: 
the entropy characteristic, 



the Alfven characteristics, 



R a' ± dx 



T dx 



R a' z du y a! y du z £^ ( a' y dB z a' z dB 



i' ± dx J~p~ \a' ± dx 



a' ± dx 



(51) 



(52) 



Ro a' z du y 
R a' ± dx 



a'y du z 



? R 



dx yfp 



a'y dB z 
a', dx 



a', dx 



(53) 



the slow mode characteristics 



£? s + = {u x + s') 



Ro 
~R 



Tdp + dT\ + ^< L du x 
p dx dx J f dx 



I'Vp 



Ro a' z dB z 
R a' ± dx 
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>, (54) 
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p dx 
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R a', dx 



a' du x 



/' dx 
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~R 
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., (55) 
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and the fast mode characteristics 



= (u x + /') 
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p dx 
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dx 
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R a' ± dx 



du x 
dx 
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., (56) 



2j- = K - /') 
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\ R a' ± dx \ R J a'± dx 
In terms of characteristics, the time derivatives (j39|) - (l45|) are then given explicitly by 
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(57) 

(58) 

(59) 
(60) 
(61) 
(62) 

(63) 
(64) 



dt Jx 

Nonreflecting boundary conditions are obtain ed by setting to ze r o the values of the characteristics for the in ward 
propagating wave s in eq uations (155)) -([M f fsee I Thompson l(| 1987ft . iThompsonl (|1990l ). iVanaiakshi et al. 1 (|1989[ ) and 
iRoe and Balsara I (|1996t) 1. 
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